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Abstract 

We propose a message-passing algorithm to compute the Hamiltonian expectation with respect 
to an appropriate class of trial wave functions for an interacting system of fermions. To this 
end, we connect the quantum expectations to average quantities in a classical system with both 
local and global interactions, which are related to the variational parameters and use the Bethe 
approximation to estimate the average energy within the replica-symmetric approximation. The 
global interactions, which are needed to obtain a good estimation of the average fermion sign, 
make the average energy a nonlocal function of the variational parameters. We use some heuristic 
minimization algorithms to find approximate ground states of the Hubbard model on random 
regular graphs and observe significant qualitative improvements with respect to the mean-field 
approximation. 
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I. INTRODUCTION 



Finding the ground-state of a quantum system can be recast as an optimization prob- 
lem by minimizing the Hamiltonian expectation over the space of trial wave functions. In 
practice, it is important for the efficiency of the variational method, to have a succinct 
representation of the trial wave functions that accurately describe the ground-state of the 



quantum system 



lJ-|3|. Nevertheless, finding the optimal variational parameters could be a 
hard computational task even in one dimension j4| since the objective function is an average 
quantity computed over the exponentially large Hilbert space of the physical system, let 
alone the complex landscape of the energy function induced by different sources of frustra- 
tions. And the problem is more serious for fermions due to the global nature of the fermion 
sign [sl. Still, the main strategy to deal with the above variational problem, is to use Monte 
Carlo (MC) method both in computing the Hamiltonian expectation and in optimizing over 
the variational parameters 6|. 

In this paper we further develop the variational quantum cavity method introduced in Ref. 
0] to study the ground-state properties of an interacting fermion system. More precisely, 
for a given instance of the variational parameters, we connect the quantum expectations to 
average quantities in a classical system of interacting particles or spins, where the interac- 
tions are related to the variational parameters. Then, instead of MC sampling, we use the 
Bethe approximation js], or cavity method in the replica-symmetric approximation jol. [lol| 
to estimate the classical expectations. Within the Bethe approximation, the probab ility 
marginals are obtained by an efficient and local message-passing (MP) algorithm jll, Il2| : 
the estimated marginals are good as long as the interaction graph is locally tree- 
classical system is in a replica-symmetric phase and it is effectively mean-fiel d [13 1 
applications of the cavity method in quantum systems can be found in Refs. \U 
may find some connections among these papers, the (statistical) dynamical mean-field theory 




20[, and density-matrix renormalization group |21 |. 

The trial wave functions can be characterized by the type of interactions included in 
the associated classical system. We usually start from the mean-field (MF) approximation 
considering only the one-body interactions, and improve on that by adding higher-order 
interactions to capture the relevant correlations. For bosons, a good estimation of the 
quantum expectations can be obtained by considering local interactions involving only a few 



number of particles ^2|. As a result, the Hamiltonian expectation is a local function of the 
variational parameters and we can again utilize the Bethe approximation to estimate the 
optimal parameters by a higher level MP algorithm 0, ■ In the case of fermions, however, 
we have to work with global interactions involving an extensive number of particles to deal 
with the global nature of the fermion sign. Consequently, the average energy becomes 
a nonlocal function of the variational parameters and we can not exploit the local MP 
algorithms to optimize over the parameters anymore. 

In this paper we take the Hubbard model and propose a class of trial wave functions with 
both local and global interactions, where the Hamiltonian expectation can be computed 
by an MP algorithm. Using some heuristic minimization algorithms we find approximate 
ground-states of the Hubbard model in random regular graphs of degree C = 3. The results 
are considerably better than the MF predictions, and are close to the exact solutions in 
small systems. For comparison we also present some results in one- and two-dimensional 
lattices. 

In the next section we give some definitions and use the mean-field approximation to 
illustrate the main points that are relevant for the following discussions. In Sec. IIIII we 
introduce the global ansatz of the wave functions and the machinery we need to deal with 
the global interactions within the Bethe approximation. The numerical data are presented 
in Sec. IIVI and finally we summarize the results in Sec. El 



II. HUBBARD MODEL IN THE MEAN-FIELD APPROXIMATION 

Consider the Hubbard model with Hamiltonian H = Hq + Hi where 

Ho = Yl ^^4t^n4iCil ^^Jt^^t + , (1) 

i i 

Hi — — ^ ^ tij{Cj^Cicr + C^^Cjfj), 
{ij)££q,a 

with index i = 1, . . . , N that labels the sites in the quantum interaction graph £g. The and 
Cicr are creation and annihilation operators for a fermion of spin a =t, i at site i. We will work 
in the occupation number representation \n) = (c|^)"'if . . . (cjy^)"^f (c^)"!-!- . . . (cjy^)"'^^|0) 
and will take the following order of the sites and spins: (1 '\,...,N t)(l l,...,N l). 
We assume that there is a path in Sq connecting 1 —)■ 2 —)■•■■—)■ representing the 
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ordering backbone. Given a trial wave function |\E'(P)) = ^^V^lz^; eP)|n) depending on a 
set of variational parameters P, we write the Hamiltonian expectation as {'^{P.)\H\'${P)) = 
EJi^in-,P)\^[Eo{n) + E,{n)]with 

Eo{n) = {n\Ho\n} = ^ 6^(71^^, n^), (2) 

i 

Ei(n)^Rei5^^1|^(n'|i/i|n)l^ J2 (-1)^— ^e.,.(n). (3) 

[ n' ^-'-^ J fo-)G^„a 

Depending on the trial wave function we obtain different expressions for eij„{n) but we 
always have ei{ni^, ra^) = Uirii^rii^ - Uiirii^ + n^). 

The goal is to consider iJ.{n;, P) = I'^/'fe P)P as a probability measure in a classical 
system and to compute the above average quantities within the Bethe approximation. The 
classical measure is, in general, represented by /i(n; P) oc Yla (paiUda) with the set of classical 
interactions Sc = {4>a{]lda)\^ = 1; • • • ) Here, da is the subset of variables that appear in 
interaction a. 

In a MF approximation, we take a factorized trial wave function including the Gutzwiller 



P) oc exp ^Kijii^riii + ^ Bi^rii^ 



(4) 



with complex parameters Ki, Bi^, and P4. This results in the following classical measure 
n{n;P) oc Hi exp [2K^ni^ni^ + Y.a'^^i^'^io-) = J^i{ni^,nij^). By superscript P, we mean 
the real part of the parameters. Given the above measure, we find 

where we defined Ana-ip = {.K*ni^ — K*nj^) + (P*^ — B*^) with f =4, and vice versa. Here 
we can easily compute the average local energies, e.g., 

2 -I— r 

{{-l)^-<k<j'^k- e,j^{n)) ^ = -tij-^-^(Pi^Pj^ + SiaSj^) Yl Wirika = 0) - Hk{nka = 1)] , 

* ^ i<k<j 

(6) 

with = e2^."+2B«+2iJ« ^ ^2Bf^ ^ ^2iJ« ^ ^ 

P.. = e^-"+^S+2BS cos(ir/ + P,i) + e^S cos(P,i), (7) 
5,. = e^-"+^S+2BS sin(ir/ + P4) + e^" sin(P4) . (8) 
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There are a few points to mention here: First, the only difference with a bosonic system is 
the sign term (— l)^i<fc<j "-fe"^. It is clear that, in the absence of this sign and for tij > 0, we 
can minimize the average energies by setting the imaginary part of the parameters to zero. 
When the sign term is present or the tij take different signs, one can show that, starting 
from real parameters, one always remains in the real subspace of the parameters following a 
gradient descent algorithm. This is true not only for the MF wave function, but also for the 
class of wave functions that we consider in this paper. Second, in the MF approximation, the 
average of the sign term is exponentially small in the number of sites k between i and j. This 
suggests that smaller average energies are obtained by maximizing the overlap between the 
ordering chosen in the trial wave function and the quantum interaction graph £g. Moreover, 
the density profile would also depend on the ordering unless the parameters are constrained 
to respect the system's translational symmetries. This artifact of the MF approximation has 
to be cured by adding interactions to the classical interaction graph £c to correlate distant 
variables along the ordering backbone. And finally, due to the sign term the average energy 
is not a local function of the parameters. This sets some restrictions on the optimization 
algorithms that we can use in minimizing the Hamiltonian expectation. 



III. BEYOND THE MEAN-FIELD APPROXIMATION: LOCAL AND GLOBAL 
INTERACTIONS 



The simplest interactions to improve the MF approximation are local two-body or Jas- 



trow interactions JijaniaUja- [25l]. It is not difficult to guess that these interactions are not 
enough to capture the sign correlations. The interaction set could, of course, be enlarged 
by adding other many-body interactions also including different types of spins. Instead, 
here, we take another approach by introducing global variables = (— l)Sfc=i '^fcT_ Then 
the sign term (— l)J^i<'=<i"'='' can be written as ^io-^io-(~l)"^'^, which is a local function of 
the global variables Q- Accordingly, we can have global one-body interactions Qia^ia 
and global two-body interactions Tijcr^ia^ja in the classical interaction graph. We call this 
set of trial wave functions the global ansatz. In general one could have interactions of type 

Notice that the above interactions do not necessarily respect the symmetries of the system. 
However, by minimizing the Hamiltonian expectation over the variational parameters, we 



get closer to the ground state of the system and, therefore, minimizing the effect of these 
asymmetries. 

In the following, we consider the global one- and two-body interactions, i.e.: 

\ i ia ia {ij)(^£c,a j 

where, for simplicity, we are going to assume Ec = £q. As a result, we obtain 

e.,.(n,i) = -UjRe {5n,.„..,oie-^"-^-^«-'^ + 5n..n,.,oie+^"-^-^«-^} , (10) 
with Ana'ip ^ given before and 

A^.7A= E Yl '^nUka^la+ Yl ^nU>^a^la- (H) 

i<k<j {kl):k<i,i<l<j {kl):i<k<j,l>j 

The average energy is computed with respect to the following classical measure: //(X) oc 
Y[i(t)i{Xi) Y[(ij)^Sc (pij{Xi,Xj) where, for brevity, we defined Xi = {n^t, ^it, ^^4; ^4-}) and 

(f)i,{X,,Xj) = 4,5^.e^-2^S-«-«-. (13) 

The indicator functions I^-^^^ check the constraints = [—1)^^'" C,{i-i)a on the global variables 
if j = z ± 1, and J^. fixes the boundary condition ^i^. = (— l)"!"^ when i = 1. To estimate the 
average energy, we resort to the Bethe approximation, writing the local marginals in terms 



of the cavity ones satisfying the following set of equations [l3|: 



fii^j{Xi) oc 4>i{Xi) Yl 5]]0ifc(Xi,Xfe)/ifc^i(Xfc) j , (14) 

where di denotes the neighborhood set of site i in £c- These are the belief propagation (BP) 
equations ll| that are solved by iteration starting from random initial cavity marginals. In 
the replica-symmetric approximation we assume there is a fixed point to the BP equations de- 
scribing the single Gibbs state of the system. The average {ei{ni^, n^))^ is simply computed 
after the local marginal fii{Xi), which is computed like fii^j{Xi) but taking all the neighbors 
into account. In the other part of the average energy, we need to compute expectations, such 
as (6<7^ja(-l)"'"5njvn,,,oi exp {-Anai^ - A^^tp))^. To get around the problem of computing 
the average of a global quantity we rewrite it as exp(F — F){$^i„Q^{—l)'^^''5nj^ni„,oie~^"'''^)fi, 



introducing the complex measure jl{X_) oc exp {—A^^ip) yu(X) and the corresponding free 
energy F. The free energy difference in the Bethe approximation is given hj F — F = 
^..(AFj — AFi) — J2{ij)(z£^{'^^ij ~ ^^ij): where AFj and AFij are the free-energy changes 
by adding site i and the interaction between sites i and j, respectively, that is, 

g-AF. ^ n ) , (15) 

Xi kedi \ Xk / 



and similarly for the complex measure [13(]. In this way, we can compute the Hamiltonian 
expectation for the above class of trial wave functions with a local message-passing algorithm 
in a time complexity of order A^^ for sparse classical and quantum interaction graphs. Note 
that a small error in estimating the classical free energies could result in a large error in 
estimating the average energy due to the exponential factor exp(F — F). 



IV. NUMERICAL SIMULATIONS 



Having the Hamiltonian expectation for an arbitrary instance of the variational param- 
eters, we need an optimization algorithm to find the optimal parameters. This is a com- 
putationally hard problem, and we have to resort to some heuristic algorithms to find an 
approximate ground state for the system. Let us start from a local minimization algorithm 
where, in each step, we fix all the parameters except in a small region of the system and 
minimize the associated energy contribution. For instance, in case Tij^j = 0, we take the 
subset {Ki, Bi^, Bi^, Qi^, 64} and minimize the following energy: 

jedi cr (kl):k<i<l o" 

The energy function is chosen as the sum of the average energies that explicitly depend 
on the subset of the parameters. The index i is selected randomly, and the corresponding 
parameters are updated. The process ends when no local update can decrease the average 
energy. 

In another algorithm, we use a population of the parameters {P_"'\a = 1, . . . ,Np} and 
update the population in order to find smaller average energies. More precisely, in each 
step, we select two sets of parameters {P_°',P!') and find the set P"^ minimizing the average 
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energy along the line AP" + (1 — A)P'' for A G [0, 1]. Then, we replace the maximal member 
of the population with P'*^ and change the position of points a and b to somewhere between 
P'^^ and the minimal member of the population P™". That is, P'* = (P"^ + p™")/2 + P" 
and P^ = (P"^ + P™") /2 + g' for some random vectors (P", P''). 

And finally, in a gradient descent algorithm, the parameters are updated as 5P = —V ^^gp'' 
for some small and positive ?7's. This means that we need to have the cavity susceptibilities 
xf-^ji^i) = ^''^^gp'^^'"' ; which cau be obtained by taking the derivative of the BP equations. 
For instance, we have 

xf^,{X^) = m^m^Si, + {XkUX^))k^^ - C,^,, (18) 

k£di\j 

if = 1, otherwise, Xi^-iji^i) = 0- Here we defined 

/ Ki _ T.Xk^ik{Xi'Xk)fJ'k^i{Xk)Xk^i{Xk) 

l^Xk Viky^ii ■^kjl^k^A^k) 

and Cf^^ is obtained by normalizatio n ^Ji'i^j{Xi)xf\.j{Xi) = 0. These are called the 



susceptibility propagation equations [28|. Similarly, we can write the equations for the 
cavity susceptibilities for the complex measure defined in the previous section. 

We can use a combination of the above algorithms to approach the optimal parameters, 
for example, the local minimization algorithm followed by the gradient descent algorithm. 
In the following, we always start from zero initial parameters or a population of parameters 
distributed randomly around zero in the case of population dynamics. The algorithms are 
repeated a number of times to get the best outcome for different realizations of the update 
process. 

Let us consider the Hubbard model on a homogeneous quantum interaction graph where 
Ui = U,Vi = u, and tij = t. In the following, we set t = 1. For simplicity, we only 
present the results obtained with real parameters; in fact by the above wave functions and 
algorithms we find the same behaviors when we take the imaginary parameters into account. 
Starting from the MF approximation, in Fig. [H we compare the average charge density in a 
random regular graph of degree C = 3 with that in one dimension; the MF approximation 
qualitatively reproduces the expected phase transitions as the chemical potential v increases 
for fixed U . For U < Uc, we observe, in turn, an empty phase, a partially filled metallic 
phase, and a completely filled insulator phase. The figure only displays densities smaller 
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FIG. 1. The average density in the mean- field approximation for the homogeneous Hubbard model 
in a random regular graph (RRG) of degree C = 3 and in one dimension (ID). The data are obtained 
by the local minimization algorithm followed by the gradient descent algorithm and repeated for 
10 different realizations of the update process. The arrows show the expected phase transition 
points for the one-dimensional system in the thermodynamic limit. The inset shows the optimal 
parameters when the wave function is given by Ki = K, Bi^ = Bi^ = B for the Hubbard model on 
the random regular graph. 

than half-filling; the other part is obtained by the particle-hole symmetry. For U > Uc, a gap 
opens up at the half-filling density, and we observe another phase transition from the metallic 
phase to an insulator phase with zero kinetic energy. Notice that, due to the exponential 
decay of the sign term in this approximation, the model on the random regular graph 
behaves nearly as the one- dimensional system. The MF approximation correctly predicts 
the empty-metal transition point where correlations are negligible and underestimates the 
metal-insulator transition point where strong correlations are responsible for the transition. 
For the parameters, we find Ki c:^ K and Bif ^ i?^ ~ i? in the metallic phase and Ki K 
and -Bjf- ~ —-84 — ±5 with a sign that changes from one site to another in the insulating 
phase. Assuming Ki = K, Bi^ = B-^, and Bi^ = B^, we can easily find the global minimum 
of the average energy for discrete parameters in a finite region of the parameter space. Up to 
the half-filling density, we find a paramagnetic solution with Bi-^ = Bi^ = B and after that, 
a ferromagnetic solution with Bi^ = — -B4 = B. Figure [1] shows how the optimal parameters 
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FIG. 2. The minimum Hamiltonian expectation computed in the MF approximation and the 
global ansatz compared with the exact ground-state energy computed by the power method, in 
one dimension (ID), and in a random regular graph (RRG) of degree C = 3. The data are obtained 
by the local minimization algorithm and population dynamics {Np = 100) followed by the gradient 
descent algorithm and repeated for 10 different realizations of the update process. 

change by the chemical potential for the paramagnetic solution. 

Going beyond the MF approximation, we always obtain smaller energies by adding global 
interactions other than the local two-body interactions. Moreover, even with the global one- 
body interactions, we obtain results that are comparable with those obtained by the global 
two-body interactions where it is more difficult to minimize the average energy, and for small 
and loopy graphs, the approximation errors cancel out the gain from the global two-body 
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Global 
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Global 
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FIG. 3. The minimum energy and average density computed in the MF approximation and the 
global ansatz in one dimension and size = 100. The data are obtained by the local minimization 
algorithm and population dynamics {Np = 100) followed by the gradient descent algorithm and 
repeated for 10 different realizations of the update process. 



interactions. Here, we present the results for the case Fjjo- = where the classical interaction 
graph has no loops, thus, the Bethe estimation of the Hamiltonian expectation is exact and 
an upper bound for the ground-state energy; this is the case even if we had the global on-site 
interactions TiC,ii^C,ii in the classical interaction graph. The global interactions in these wave 
functions can be considered as effective ones representing higher-order global interactions. 
The associated wave functions are expected to describe the physical state of the system for 
U ^ t well. In Fig. [2] we compare the minimum energies with the exact solutions on small 
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FIG. 4. The minimum energy and average density computed in the MF approximation and the 
global ansatz in a random regular graph of degree C = 3 and size N = 500. The data are obtained 
by the local minimization algorithm and population dynamics {Np = 100) followed by the gradient 
descent algorithm and repeated for 10 different realizations of the update process. 

systems computed by the power method using an infinitesimal imaginary time evolution 
1 — Ht to minimize the Hamiltonian expectation 29|. Indeed, the power method can also 
be implemented within the variational formalism where, in each step, one has to project the 
change in the wave function onto the space of the variational parameters; see Ref. jsol for 
more advanced methods. 

Figures [3] and H] display the average charge density for some larger system sizes. In one 
dimension we are very close to the expected phase transition points in the thermodynamic 
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FIG. 5. The spin and charge density profiles in the global ansatz for different chemical potentials 
V — U /2 = —0.2, —0.5, —1 (from top to bottom) and [/ = 5 in a random regular graph of degree 
C = 3 and size = 100. The solution is obtained by the population dynamics [Np = 100) followed 
by the gradient descent algorithm and repeated for 10 different realizations of the update process. 

limit. And in random regular graphs, we observe a shift to smaller chemical potentials for 
the empty-metal transition point and a shift to larger chemical potentials for the metal- 
insulator transition point, with respect to the one- dimensional model. One can attribute 
these shifts to the larger connectivity (here, C = 3) of the random graphs that provide more 
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FIG. 6. The average energy and charge density obtained in the MF approximation on the 2D square 
lattice. The data for paramagnet to ferromagnet (P-F) and antiferromagnet (P-AF) transitions 
are obtained by searching for the optimal parameters {Ki = K,Bi^ = B^,Bi^ = B^) and {Ki = 
K, Bi^ = Bj^ = B, Bi^ = Bj^ = B) (for i in odd and j in even sub-lattices), respectively. 

degrees of freedom to the particles. 

To obtain a picture of the wave functions, in Fig. [5|, we show the spin and charge density 
profiles along the ordering backbone of our representation. In the insulator phase, the spin 
densities are not frozen on or 1 anymore as happens in the MF approximation. In addition, 
we observe spin density oscillations that were absent in the MF approximation and nonzero 
kinetic energies in the insulating phase. Before the metal-insulator transition, we observe 
some charge density holes separating different antiferromagnetic regions. The global one- 
body parameters Gjo- alternate between positive and negative signs and are different for the 
spins up and down. 

Finally, we report some preliminary results for the Hubbard model on the 2D square 
lattice. Figure [6] shows the average energy and density obtained by the MF approximation 
when we allow for a phase transition from the paramagnetic (P) phase to the ferromagnetic 
(F) and antiferromagnetic (AF) phases. The former transition happens before the latter 
one, and the difference between the two energies decreases as one approaches the half-filling 
density where the AF phase has a smaller energy. Adding the global one-body interactions, 
we find a transition from the P phase to a mixed phase of F and AF regions, see Fig. O the 
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FIG. 7. The spin density profiles in the global ansatz for chemical potentials v — [7/2 = —5, —10 
(from top to bottom) and ?7 = 20 in ID square lattice of size = 10 x 10. The solution is obtained 
by the population dynamics [N^ = 100) followed by the gradient descent algorithm and repeated 
for 10 different realizations of the update process. 



system is more ferromagnetic (antiferromagnetic) for smaller (larger) chemical potentials. 
Figure [H] displays the charge and spin profiles close to the half-filling density. We observe 



tendencies towards the holes condensation 



3l| where a ferromagnetic phase of small density 



is separated from an antiferromagnetic phase of higher density. 
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In Fig. IHl we display the time te we need to compute the Hamiltonian expectation and the 
local average quantities in the 2D square lattice with the global one-body interactions, given 
the set of parameters Ki, Bi^, and Gjcr. The computation time for updating the population 
of the parameters (one sweep) is Nptg, and, in practice, we need a few hundred sweeps of the 
updates to reach the stationary state. We are reminded that, as long as one only considers 
the global one-body interactions, the algorithm gives the exact Hamiltonian expectation 
and, therefore, a good upper bound for the ground-state energy. For comparison, we also 
show the computation time for a diluted lattice where a small fraction of the global two-body 
interactions are present in addition to the ones in the ordering backbone. 



V. DISCUSSION 



The Hamiltonian expectation can be computed by an efficient and distributive message- 
passing algorithm, which is asymptotically exact on random and sparse interaction graphs. 
To obtain a good estimation of the average fermion sign, we have to work with the global 
interactions in the classical system. Unfortunately, this makes the average energy a non- 
local function of the variational parameters, resulting in an optimization problem that is 
not amenable anymore to local message-passing algorithms. Moreover, the performance 
of any optimization algorithm strongly depends on the quality of the estimated average 
energy or the approximation method that is used to take the average of the energy func- 
tion. The study can systematically be improved in both directions b y c onsidering more 



accurate inference algorithms using generalized Bethe approximations [32|, |33| or incorpo- 
rating replica-symmetry-breaking and more sophisticated optimization algorithms to find 
the optimal variational parameters. 
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